c 



Temperature and polarization CMB maps from primordial non-Gaussianities of the 

local type 

Michclc Liguori , Amit Yadav 2 , Frode K. Hansen 3 , Eiichiro Komatsu 4 , Sabino Matarrese 5 , Benjamin Wandclt 2 

1 Department of Applied Mathematics and Theoretical Physics, 
Centre for Mathematical Sciences, University of Cambridge, 
Wilberfoce Road, Cambridge, CB3 OWA, United Kingdom 
2 Department of Astronomy, University of Illinois at Urbana- Champaign, 1002 W. Green Street, Urbana, IL 61801 
3 Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, 0315 Oslo, Norway 
4 Department of Astronomy, University of Texas at Austin, 2511 Speedway, RLM 15.30 6, TX 78712 and 
° Dipartimento di Fisica "G. Galilei", Universitd di Padova and INFN, 
Sezione di Padova, via Marzolo 8, 1-35131, Padova, Italy 
(Dated: February 1, 2008) 

The forthcoming Planck experiment will provide high sensitivity polarization measurements that 
will allow us to further tighten the /nl bounds from the temperature data. Monte Carlo simulations 
of non-Gaussian CMB maps have been used as a fundamental tool to characterize non-Gaussian 
signatures in the data, as they allow us to calibrate any statistical estimators and understand the 
effect of systematics, foregrounds and other contaminants. We describe an algorithm to generate 
high-angular resolution simulations of non-Gaussian CMB maps in temperature and polarization. 
We consider non-Gaussianities of the local type, for which the level of non-Gaussianity is defined 
by the dimensionless parameter, /nl. We then apply the temperature and polarization fast cubic 
statistics recently developed by Yadav et al. to a set of non-Gaussian temperature and polarization 
simulations. We compare our results to theoretical expectations based on a Fisher matrix analysis, 
test the unbiasedness of the estimator, and study the dependence of the error bars on /nl. All our 
results are in very good agreement with theoretical predictions, thus confirming the reliability of 
both the simulation algorithm and the fast cubic temperature and polarization estimator. 

I. INTRODUCTION 

Small, but non-vanishing non-Gaussianity of primordial cosmological perturbations is a general prediction of in- 
flation. The amplitude of the expected non-Gaussian signal is model-dependent and can vary by many orders of 
magnitude from one inflationary scenario to another. For example, the non-Gaussian signatures produced by single- 
field slow-roll inflation models are tiny and far below the present and forthcoming experimental sensitivity [l|, [2| . On 
the other hand many other scenarios predict a level of non-Gaussianity that is within reach of present and forthcoming 
experiments like WMAP and Planck (see e.g. d i, H, i, 0, H, S For this 

reason an experimental detection of 

non-Gaussianity would rule out the simplest scenarios of slow-roll inflation. More in general, experimental bounds on 
primordial non-Gaussianity allow us to significantly constrain different scenarios for the generation of perturbations 
in the context of primordial inflation. 

Primordial non-Gaussianity from inflation can be described in terms of the 3-point correlation function of the 
curvature perturbations, $(k), in Fourier space: 

(*(ki)*(k a )$(ks)> = (27r) 3 S^(k 1 +k 2 + k 3 )F(k 1 ,k 2 ,k 3 ) . (1) 

Note that <£> is the curvature perturbation during the matter era, and temperature anisotropy in the Sachs- Wolfe 
limit is given by AT/T = — $/3. Depending on the shape of the function F(k±, k 2 , k 3 ), we can divide non-Gaussianity 
from inflation into two different classes: local non-Gaussianity, where F is large for squeezed configurations (i.e. 
configurations in which k\ << £2,^3), and non-local non-Gaussianity of the equilateral type, where the largest 
contributions come from modes with k\ ~ k 2 ~ &3- The former kind of non-Gaussianity can be produced in models 
where primordial perturbations are not generated by inflaton itself but by a second light scalar field (like e.g. in 
the curvaton model). The latter comes from single field models with a non- minimal Lagrangian containing higher 
derivative operators. In this paper we will focus on non-Gaussianity of the local type, where the primordial curvature 
perturbation $ can be described in terms of the following real space parameterization: 



$(x) = *i(x) + /nl ($i(x) - (*|(x)» 



(2) 
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In the last formula /nl is a parameter that defines the amplitude of the primordial non-Gaussian signal. Our 
previous statement about the detectability of non-Gaussian signatures from inflation can be precisely quantified in 
terms of this parameter. In standard scenarios of single-field slow roll inflation /nl is generally predicted to be very 
small and undetectable (~ 10 -2 at the end of inflation, ~ 1 when second order perturbation theory after inflation is 
taken into account) whereas other scenarios, like the curvaton or variable decay width models, can naturally give rise 
to relatively large values of /nl (/nl ~ 10). This justifies the claim that an experimental detection of /nl would rule 
out the simplest single-field inflationary paradigm and allow us to put significant constraints on the other inflationary 
scenarios. 

The best way to put experimental bounds on /nl is to look for non-Gaussianities in CMB anisotropics (but it 
has been recently pointed out that future deep galaxy-surveys and 21 cm background measurements could provide 
promising results [IH [3, EH ) . The most stringent constraints on /nl so far come from measurements of the CMB 



15, [16|]. This constraint 
13 showed that WMAP 



angular bispectrum on the WMAP temperature data —36 < /nl < 100 (95% ci.) \l£ 
corresponds to a 1 — a error of A/nl = 34. A Fisher matrix analysis by the authors of 
will in principle be able to reach A/nl = 20, while the forthcoming Planck satellite can achieve A/nl = 5. This 
means that Planck will be sensitive to the level of non-Gaussianity predicted by a vast range of different inflationary 
models. We can improve this constraint further by including the polarization data. For WMAP all the non-Gaussian 
information is basically contained in the temperature data, due to large errors in polarization measurements. Planck, 
on the other hand, will characterize polarization fluctuations with high accuracy. This will allow us to exploit the 
additional information contained in polarization data and to gain a further factor of order 2 in A/nl, thus yielding 
A/nl — 3 18]. A crucial step in order to exploit all the information contained in the future Planck dataset is then to 
extend the tools previously developed for temperature non- Gaussianity in order to include polarization. This program 
has been recently started by the authors of [13], where the fast cubic statistic used to analyze WMAP temperature 
data [20l . l2l| was taken as a starting point to build an optimal cubic estimator that is sensitive to a combination of 
temperature and polarization primordial fluctuations. In this paper we will extend the non-Gaussian analysis toolkit 
in order to include the second fundamental element: Monte Carlo simulations of primordial non-Gaussian polarized 
CMB maps. 

In section [TT] we will summarize the original algorithm and describe its extension to polarization. We will then 
apply the fast cubic statistic of [191 ] to a set of polarized non-Gaussian maps. In this paper we will manly focus our 
attention on map generation, so the purpose for applying the estimator is mainly to check the reliability of the final 
maps. This will be done by comparing the final outputs to theoretical predictions in ideal conditions. However in a 
forthcoming publication we will describe how we actually used the maps in order to test, calibrate and optimize the 
estimator. 



II. GENERATION OF POLARIZED NON-GAUSSIAN CMB MAPS 



Realistic simulations of non-Gaussian CMB maps are indispensable tools for measurements of non-Gaussian signals 
in the data, as they allow us to test and calibrate estimators and also to include and study all the spurious non- 
Gaussian signals introduced by contaminants like foregrounds, secondary anisotropies, instrumental noise and so 
on. 

The first simulations of temperature maps with primordial non-Gaussianity from /nl were carried out by Komatsu 
et al., and used extensively to study Gaussianity of the WMAP data [IH as well as non-trivial topology of the universe 
[22I ]. Then, Liguori et al.[23[ have succeeded in increasing the computational speed, reducing the memory requirement 
and, most importantly, improving accuracy of the simulated temperature maps. We take this new algorithm developed 
in [23j as a starting point. 

Our starting point is the relation between the primordial curvature perturbation $ and the CMB multipoles af m 
via radiation transfer functions A^. 

/d 3 k 
<f,(k)Yt m (k)A? (k) , (3) 

where X refers to either the temperature component T or the polarization component E. 

The kind of non-Gaussianity we are considering has a very simple form in real space, where it is local and the 
non-Gaussian part of the curvature perturbation is simply the square of the Gaussian part (see formula [2|) . For this 
reason it is convenient to work in real space and define the real space transfer functions Af (r) as: 

Af (r) ee - / dkk 2 Jt (kr)Af (k) , (4) 
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FIG. 1: CMB angular power spectra extracted from 10 simulations (triangles) are compared to the theoretical ones computed 
with CMBfast for the same model (solid black lines) . The cosmological parameters are = 0.042, f2 c dm = 0.239, Qr, = 0.719, 
h — 0.73 n — 1, and r = 0.09 (same for all the following figures, unless otherwise stated). 



where ji(kr) is the spherical Bcsscl function of order I. It can be shown that A^(r) links the primordial curvature 



perturbation $(r) in real space to the af m through the following relation [14l, |2 



x 

a lm 



drr 2 Af(r)<S> lm (r) 



(5) 



In this last formula we have introduced the quantities < I ) £ m (r), which represent the spherical harmonic expansion 
multipoles of the curvature perturbation $(r, r) on a shell of given radius r. In formulae: 



$*m(r) - / dSlfY tm (fMr,f) . (6) 

We define the radius r as r = c(ro — r), where c is the speed of light and tq — r is the lookback conformal time. 
The radius r varies from the origin r — to the present time cosmic horizon r — ctq. The radii in which $^ m (r) must 
be generated depend on the features of the real space transfer function Af (r) in equation (JSJ) . We will come back to 
this shortly. 
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FIG. 2: Temperature (bottom panel) and polarization (upper panel) transfer functions at high i at last scattering. 



Let us assume for the moment that we have been able to numerically generate the Gaussian part of the curvature 
perturbation multipoles Qf m (r) for the chosen set of radii. Starting from here we can now generate the non-Gaussian 
part QfJ^ (r) in the following way. First of all we harmonic transform (r) to get the gaussian part of the curvature 
perturbation in real space: 



$ L (r,f) = ^^^ m (r)y, ro (f). 



(7) 



l m 



Then we square <&h(r, r) to get the non Gaussian part of the curvature perturbation on each sampled spherical shell: 
^NLfj, f) = $l0">^) — (^t( r ^))- We then calculate the multipoles of this non-Gaussian part through a backward 
harmonic transform: 



$^(r) = / dn^ Nh (r,r)Y em (f) 



(8) 



Having computed Qf m (r) and ^^(r) we can finally obtain the Gaussian and non-Gaussian part of the CMB 
multipoles, af^ and af^ h respectively, by applying formula ([5]): 



X,NL 



drr 2 A?(r)^ m (r) 
drr 2 Af(r)$Z(r) 



(9) 
(10) 
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FIG. 3: Temperature (bottom panel) and polarization (upper panel) transfer functions at low i (reionization and late ISW 
contributions are visible). The oscillations visible in the plots are little numerical artifacts which have negligible impact on 
the final results. We have explicitly checked this by increasing the resolution in the k and r-grid by factors of 2 and 4 without 
noticing any improvement in the accuracy of the final Ci, that can be already reconstructed well using the sampling chosen in 
the paper (see fig. [1} 



A CMB map for a chosen value of /nl can then be obtained simply by summing af^ + /NLa^ NL - This means 
that with a single generation of af m and a^ NL it is possible to generate maps for any value of /nl- 

We are still left with one problem unsolved i.e. how do we generate the Gaussian curvature perturbation multipoles 
&f m (r) ? This issue is complicated by the fact that curvature perturbation multipoles are correlated in real space. 
The obvious solution would be to generate curvature perturbations in Fourier space, 3 > (k), Fourier transform back 
to real space to obtain <&(x), change the coordinates from Cartesian to polar to obtain $(r, h), and finally harmonic 
transform to obtain <&£ m (r). This is the original approach taken by [14| . which is computationally quite expensive. 
Also, the coordinate transformation from Cartesian to polar limits accuracy of the maps, especially at high multipoles. 

A novel approach developed in [23[ solves this issue by generating $^ m (r) directly, without ever worrying about the 
coordinate transformation. It has been shown in [23| that the <&i m {k) and ^ m (r) are related by a spherical Bessel 
transform: 



*/m(r) = J dkk 2 Je {kr)<S> lm {k) . (11) 

The problem with this expression is that the Bessel functions oscillate very rapidly. This implies that, for each 
(£,r), the integral above must be sampled in many different k in order to attain sufficient accuracy, thus making the 
computational cost of such an algorithm prohibitive. A much more convenient solution was found in [23|; the idea is 
to start with a set of Gaussian independent "white noise" coefficients ri£ m (r) characterized by the following correlation 
function: 
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FIG. 4: Temperature transfer functions at high i and r corresponding to the epoch of reionization. Polarization transfer 
functions at large £ are zero in this range. 



(ne imi (ri)n* e27n2 (r 2 )) 



} (n-r 2 ) 



°e 1 °m 1 i 



(12) 



it can be now shown that Gaussian curvature perturbation multipoles (fc^ (r) with the right correlation properties 



can be obtained through a convolution of the ne m coefficients with suitable "filters" Wf. 



(13) 



where the functions We are defined as 



W e (r,n) = - dkk 2 VMk)je{kr)h(kn) , 

7T 



(14) 



and P®(k) is the power spectrum of the primordial curvature perturbation <&L(k). As depicted in Fig. [6l [7] the 
filter functions We are smooth. Moreover, as also suggested by the Limber approximation applied to equation (|14p . 
We(r, ri) is narrowly peaked around r when I > 10. This allows to sample the integral (|13p in much less points than 
it would be required for the Bessel transform (jlip . thus making the problem computationally feasible. Obviously the 
problem of sampling a highly oscillatory integrand has not disappeared completely but it has been reduced to the 
generation of Wi(r, ri). A trick here is that the filters Wi(r, ri) can be pre- computed and stored once and for all for a 
given cosmological model and their calculation does not enter in the actual Monte Carlo simulation algorithm. The 
same argument applies to the radiation transfer functions Ag(r) defined in (|U). 
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Region 


Bounds 


Ar 


N. of shells 


Recombination 


12632 Mpc < r < 13682 Mpc 


3.5 Mpc 


300 


Reionization 1 


10007 Mpc < r < 12632 Mpc 


105Mpc 


25 


Reionization 2 


9377 Mpc < r < 10007 Mpc 


35Mpc 


18 


Low redshifts 


Mpc < r < 9377 Mpc 


105Mpc 


89 



TABLE I: Sampling of the r-coordinate in different regions of the simulation box. Different intervals must be sampled with 
different resolutions, according to the radiative transfer physics described in section III Al 

Non-Gaussian temperature maps produced with the algorithms described in this section had been already described 
in [23| . Adding polarization to those maps is conceptually straightforward: all one needs to do is to replace X = T 
with X = E in the previous expressions. This amounts to generating the primordial curvature perturbation in exactly 
the same way for temperature and polarization maps and finally to use polarization transfer functions in place of 
temperature transfer functions in the line of sight integral ([5]) in order to get af m . Despite its conceptual immedi- 
ateness, including polarization in the maps is not technically straightforward. The reason is that CMB polarization 
is produced by different physical mechanisms with respect to those producing CMB temperature anisotropics. The 
polarization transfer functions Af (r) present then several differences with respect to Aj(r) and must be sampled in a 
different way, thus changing sampling regions and discretization of the r-coordinate which appears in <&i m (r), A^(r), 
Wi{r, ri). These technical details will be illustrated in the following two sections. 

A. Real space transfer functions 

The cosmological model we chose to generate our non-Gaussian maps is characterized by the following parameters: 
ftcdm = 0.239, fib = 0.042, = 0.719, r = 0.09, h — 0.73. We considered both a scale invariant primordial 
spectral index n — 1 and n — 0.95, the latest one being the WMAP 3-years best-fit value [l5j]. Starting from 
these parameters we generate and extract the Fourier space radiation transfer functions Af (k) from a Boltzmann 
integrator, like for example CMBfast, and then make the integral ([4]) to get Af{r). The behavior of Af (r) reflects the 
underlying temperature and polarization CMB physics. In Fig. [2j we plot the real space temperature and polarization 
transfer functions for several different values of I > 20. For the model under examination the conformal time at last 
scattering, defined as the peak of the visibility function, is r* ~ 277 Mpc (c = 1) while the present cosmic horizon is 
r ~ 13682 Mpc. We thus expect most of the signal to be generated at r* = t — t* ~ 13400 Mpc, consistently with 
what shown in the figure. Despite being smaller, contributions at lower redshifts cannot be neglected. We know that 
both reionization and the late integrated Sachs Wolfe effect produce significant contributions, especially at low £'s. 
The reionization signal is particularly important for polarization, as it produces the observed bump at low £'s in the 
polarization spectrum. This is reflected in the behavior of the temperature and polarization transfer functions at low 
I in the post-recombination region, accordingly to what depicted in Fig. 2] and Fig. [31 According to the radiative 
transfer physics contained in Ae(r), the last scattering surface must be sampled using a large number of points in 
order to accurately reproduce the acoustic oscillations in the CMB spectrum, while in the low redshift region a good 
accuracy can be reached with a coarser sampling. More details about the sampled regions and intervals are in table 
HI the idea was to refine the r-grid until a good accuracy in the final CV from the simulated map was reached (see Fig. 
[T]). However further sampling optimization in order to improve the speed of the algorithm is probably possible; an 
algorithm aimed at this kind of optimization is described in [2~ij in the context of bispectrum estimation. 



B. Filter functions 

After generating the real radiation transfer functions and fixing the radial coordinate grid, the We(r,n) functions 
defined in (jT5J) must be generated for each value of £, r. Due to the highly oscillatory nature of the Bessel functions 
appearing in the definition of Wi(r,r\), a large number of points is required when sampling the integrand. This 
makes the numerical computation of Wg(r,ri) quite slow. However, as we were already stressing above, this is not 
a problem as the We (r, r\ ) functions are pre-computed and stored before the actual Monte Carlo map generation. 
When computing Wi it is useful to make the simple substitution t = kr in the integrand of fl!4[) . This substitution 
yields: 
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W i [r,r 1 )=2-Kr = ^h(^) , (15) 

where we have defined: 



h{x) = J dtt^ ji(t)ji{tx) . (16) 

From the last formulae we see that We(r, ri) actually depends only on the ratio r\/r and not on r\ and r separately. 
This allows to reduce the dimensionality of the problem and thus to speed up the calculations. 

In figures [5] and [7] we plot some Wg functions for different values off, r, r\. As expected, W^(r, r x ) approximates 
a Dirac delta function centered on r with increasing accuracy for larger and larger £. So for I > 10 the coordinate 
ri needs to be sampled in a narrow region centered around r. On the other hand, for low values of £, Wi(r,ri) is 
non-negligible over a broad range. Thus a coarser sampling over a larger r-± interval is required in this case. To 
check the accuracy of the numerical computation of Wg(r, r{) it is useful to compute the angular power spectrum of 
$£ m (r) on a given spherical shell. Starting from formula (|13p . and using the correlation properties of the coefficients 
ri£ m (ri) described by eqn. IT2|) one gets: 

(^m.W^kW) = ^S ei i 2 S mi , n2 J dr 1 dr 2 rlrl[{n iimi {r 1 )n* l2m2 {r 2 )) X 

x W tl (x,n)W (2 (y, r 2 )] 

J '— lir 

2 



dr x dr 2 r 1 r 2 3 W *i ( x > r i) w t2 {V^i) 

= -S ei i 2 5 mi m 2 / drirlW ei (x,ri)We 2 (y,r 2 ) , (17) 



which immediately yields: 



(l$LMI 2 ) = f / dnriwiir,^) . (18) 

Alternatively it is possible to use the following formula for the correlation function [23j |: 

(x)®f; m2 (V)) = ~%S% J dkk 2 P*(k)j ei {kx)j i2 (ky) , (19) 

to find: 

(\®L(r)\ 2 ) = lj dkk 2 P(k)f t (kr) . (20) 

For a primordial curvature perturbation power spectrum described by a power law expression, P(k) — Ak n ~ 4 , and 
using a well-known formula for the Sachs- Wolfe effect, one finally gets: 

„ r , Xl9l 2 n - 3 Ar 1 - n TU+ % - ±) T (3 - n) / s 

For a scale invariant primordial power spectrum one obtains, as expected, Id^^MI 2 )! ^ + 1). As we were 
anticipating above, one can use formulae (fT5| and (j2~lj) to test the $£ m (r) power spectrum on different shells and the 
normalization of Wg (r, r\ ) . Results from our simulations are shown in picture [5j 
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Noise 


Sky-cut 


(f NL ) 


^maps 


^fisher 


No 


No 


102.5 


11.1 


6.9 


Homogeneous 


No 


104.5 


15.8 


11 


Homogeneous 


fsky = 80 


105.2 


25.7 


12.2 



TABLE II: Results obtained from the application of the fast temperature + polarization cubic statistics of [TS| to a set of 300 
non-Gaussian maps with an input /nl of 100. First column describes the noise properties of the map, second column is the 
adopted sky-cut, third column is the average /nl measured by the estimators, fourth column is the measured /nl standard 
deviation, fifth column is the expected standard deviation from a Fisher matrix analysis (i.e. neglecting corrections from the 
non-Gaussian part of the multipoles). 
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FIG. 5: Angular power spectrum of the Gaussian curvature perturbation multipoles <&£ TO (r) obtained by averaging over all the 
spherical shells of a given simulation. In this example we consider a spectral index n = 0.95 and divide |$£ m (r)l by V r( 1_n ) 
in order to make the normalization of the spectrum independent of the shell radius before averaging. We compare the results 
extracted from our simulations (red triangles) to the expected shell power spectrum obtained from formula (I21[l . (blue line) 



III. FAST CUBIC STATISTICS AND NON-GAUSSIAN MAPS 

In order to test our algorithm we applied the temperature + polarization fast-cubic statistics described in (l9| 
to a set of 300 non-Gaussian simulations obtained from the cosmological parameters ilt = 0.042, fl c d m = 0.239, 
fli = 0.719, h = 0.73 n = 1, r = 0.09. In figure [S] we show a temperature and a polarization intensity map extracted 
from this set. 

When skycut is included a non-trivial correlation between large and small £ is introduced. This correlation in turn 
produces a leakage of power from high to low multipoles which tends to bias the estimator. This effect has been 
accurately studied in [191 ] , where it has also been shown that removing the lowest multipoles from the analysis allows 
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FIG. 6: Filter functions Wt{r, n) plotted as a function of ri for two different fixed values of r. Here we consider low 1-values 
I < 10, for which the Wi(r, ri) are different from zero and must therefore be sampled in a large n region. At high I these 
functions become more and more peaked around r, as shown in Fig. [7] 



to circumvent this problem without a significant loss of signal. For this reason the first 30 multipoles were not used in 
our analysis when a skycut was considered. The exact £ m in was determined by preliminary applying the estimator to 
a set of Gaussian simulation and estimating its variance as a function of £ m in- We considered different sky cut levels 
and accounted for the presence of homogeneous noise. Our results are summarized in table HT1 . 

Our computation provides evidence for the unbiasedness of the estimator but shows at the same time a discrepancy 
between the calculated error bars and Fisher matrix based expectations (note that these expectations are obtained 
at zeroth order in /nl, thus neglecting /NL-dependent terms in the three point function). These discrepancies are 
/nl dependent: for small undetectable /nl we find a good agreement between Fisher matrix estimates and our 
results whereas increasing values of /nl produce larger and larger differences. This effect had been predicted and 
explained by the authors of [l6|. It arises from /nl dependent correction terms in the variance of the estimator. 
These terms become important when /nl is detected at several sigma. The comparison of our results with those in 
[l6| is necessarily approximate because the latter were obtained in the flat-sky approximation and ignoring radiation 
transfer functions. However we can still cross-check for a qualitative agreement between the two results. Using the 
above approximations, the /NL-dependent formula describing the estimator variance is: 



= {a ) /NL=0 1 + — , (22) 



7r In Npix 



where (a 2 ) / NL =o is the estimator variance in the Gaussian case (i.e. the variance estimated from the Fisher matrix), 
A is the amplitude of primordial perturbations and N p i X is the number of pixels in the map. To simplify the notation 
we define ctq = (<7 2 )/ NL= o- Following (l6l | we consider an /nl detection at na®. From the formula above: 
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FIG. 7: Filter functions We(r, ri) plotted as a function of ri for two different fixed values of r. As I gets larger, the We(r, n) 
becomes more and more narrowly peaked around r. 



2n 2 al 
7T In 2 N n , 



(23) 



We then find the expected relative correction to the variance as: 



- 1 = 



2n 2 



tt In 2 7V„ 



(24) 



In our analysis we have N p i X — 3145728 (HEALPix nside = 512) and Co = 6.9 for the case without sky-cut or noise 
(see second line of table [TTJ) . We have an input /nl of 100, so this corresponds to n — 14.5. Plugging this numbers 
into the equation above we obtain a relative correction of 0.6 which is about one third of the observed a 1 /<j$ — 1 = 1.6 
but in qualitative agreement considering the approximations contained in eqn. (|24p. For large enough /nl eqn. 



also predicts the variance of the estimator to decrease as 1/ln N p i X ~ 1/ \nl max , much slower than the Fisher matrix 
forecast of a ~ l/£ m ax- We explicitly tested this prediction on sets of simulated maps with different /nl , N p i X , i rna x 
and we found a good agreement between theory and simulations, as depicted in figure [9] Thus the results obtained 
analytically in [lfj under several simplified assumptions are confirmed by our numerical approach, which works in 
full-sky and includes radiation transfer functions. 
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FIG. 8: Left column: temperature an d polariza tion intensity Gaussian CMB simulations obtained from our algorithm. Po- 
larization intensity is defined as I = a/Q 2 + U 2 where Q and U are the Stokes parameters. Right column: temperature and 
polarization non-Gaussian maps with the same Gaussian seed as in the left column and /nl = 3000. The reason for the choice 
of such a large /nl is that we wanted to make non-Gaussian effects visible by eye in the figures. The cosmological model 
adopted for this plots is characterized by: £lb ~ 0.042, £l c dm = 0.239, Ql = 0.7f9, h — 0.73, n = 1, r = 0.09. Temperatures are 
in mK. 



IV. COMPUTATIONAL REQUIREMENTS AND POSSIBLE APPLICATIONS 

Our algorithm takes about 3 hours on a normal PC to generate a map with £ max = 500, N p i X ~ 10 6 , corresponding 
to an analysis at WMAP angular resolution. The most time consuming part is the computation of the harmonic 
transforms required to generate <I>^ from $^ TO . As we generate the primordial curvature perturbation in about 400 
spherical shells we need to make 400 calls to the HEALpix synfast and anafast subroutines respectively. Thus we can 
roughly quantify the CPU time for a non-Gaussian simulation at a given resolution as the time required to produce 800 
Gaussian maps at the same resolution. It is thus clear that the generation of maps at the resolution achieved by Planck 
constitutes a very intensive computational task and requires a parallelization of the algorithm. Only the temperature 
version of the code has been parallelized so far, enabling us to generate a map at £ max = 3000, nside = 2048 in 
about 2 hours on 60 processors. A set of 300 temperature maps with this angular resolution has been generated 
and tested. Extending the parallel code in order to include polarization should be straightforward, because all the 
sampling-related problems have been already solved for the serial version of the algorithm presented in this paper and 
including polarization transfer functions is trivial. The total CPU time to generate a map is going to be unchanged 
with respect to the temperature-only version, because the primordial curvature perturbation generation scheme is 
identical and the total number of shells is basically the same. We would like to note here that a different algorithm 
has been proposed for the generation of non-Gaussian maps in [24j . This algorithm can generate maps with a given 
two and three point function but does not reproduce the higher order correlation functions predicted by the model. 
By making this approximation, the authors of [24| are able to dramatically speed up the computation (~ 3 minutes 
for a map at £ m ax — 1000 on a single processor). In the limit of weak non-Gaussianity citenotel neglecting higher 
order correlation functions should be a good approximation. In particular it has been explicitly shown in [l6| that no 
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FIG. 9: Error bars estimated from different sets of simulations including various £i ma x and input /nl. The error bars are 
compared to the corresponding Fisher matrix forecast. As explained in the text, an /NL-dependent correction to the estimator 
variance make the error bars to scale as l/\n£ ma x instead of l/£ ma x when /nl is large enough to produce a several sigma 
detection at a given angular resolution. 



additional information on /nl can be added by applying estimators based on higher order correlators. This conclusion 
is strictly related to the presence of /NL-dependent correction terms in the variance of the local bispectrum. Note 
however that these terms have originally been studied in flat-sky approximation and neglecting transfer functions. As 
an application of our algorithm, in the previous section of this paper we have explicitly cross-checked the results of 
16] using our simulations which are full-sky and account for radiative transfer (31| 

We would also like to stress that being able to correctly reproduce higher order correlation functions in the sim- 
ulations was fundamental in order to make this test. The reason is that what we are studying here is actually an 
/NL-dependent correction to the 6-point function (bispectrum variance) coming from a product of the 2-point function 
with the 4-point function (see again [l6| for further details). 

Another obvious application for these simulated maps is given by the possibility to use them in order to test and 
calibrate not only the bispectrum but any kind of estimator (like e.g. Minkowski functionals, wavelets and so on). In 
particular the analytical formulae of the Minkowski functionals recently derived by [25[ may be compared with our 
simulations of the temperature maps. Our preliminary investigation shows a very good agreement, which gives us 
further confidence in the accuracy of the simulated temperature maps. Despite the optimality of the bispectrum just 
discussed above, using different estimators is still important, especially in view of a possible /nl detection by Planck. 
Alternative estimators should in fact be used in this case in order to cross-validate such detection. 

Furthermore, it is interesting to notice that the algorithm we are describing is not only able to generate non-Gaussian 
CMB maps, but it also produces maps of the primordial curvature perturbation $(r, f), sampled in the relevant radii 
for the generation of the final CMB signal. This allows us to apply and test tomographic reconstruction techniques 
of the curvature perturbation like those proposed in (26[. This will be the object of a forthcoming publication [27| . 
Finally we would like to observe that the same elegant r-sampling optimization technique introduced in [24j can 
be implemented in our case in order to drastically reduce the number of radii in which the primordial curvature 
perturbation must be evaluated. Following the results of dU, a good accuracy in the final maps should be obtained 
using only 20 spherical shells in our code after this optimization. As we are now using 400 shells, we estimate a speed 
improvement of a factor ~ 20. In this way the parallel version of the algorithm should allow the generation of a map 
at full Planck resolution in ~ 10 minutes against the present 2 hours. For this reason CPU time does not seem to be 
a problem and tests of non-Gaussianity at Planck angular resolution using our algorithm are perfectly feasible. 
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V. CONCLUSIONS 

In this paper the algorithm for the generation of non-Gaussian primordial CMB maps originally introduced in (23| 
has been generalized by including a polarization component in the simulations. Using this generalized algorithm we 
have produced a set of 300 temperature and polarization maps at WMAP angular resolution. We have then analyzed 
these simulations using the fast cubic temperature + polarization statistics recently introduced by the authors of 19]. 
We have verified that we can extract the correct input /nl from the maps, thus checking at the same time both the 
unbiasedness of the estimator and the reliability of the simulations. We also studied the estimator variance on different 
sets of maps including various angular resolutions and input /nl- We found that an /NL-dependent correction to the 
estimator variance induces a discrepancy between the error bars extracted from the simulations and the Fisher matrix 
estimate of the same error bars at /nl = 0. We therefore confirmed previous findings by the authors [l6j]. At the 
same time, differently from previous approaches, our numerical Monte Carlo analysis allowed us to work in full sky 
and account for radiation transfer functions. We finally discussed future applications of our simulations, which will 
include a detailed analysis of non-Gaussian temperature and polarization simulations at Planck angular resolution. 
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This limit is verified in our case as non-Gaussianity from inflation is small. 

As a subject of future work, all these checks will be repeated by taking into account possible contaminant effects, like e.g. 
foreground residuals, second order anisotropies, systematics, map-making effects and so on, in order to study their impact 
on the estimator. 



